R Notebook

Libraries

Load data

I am using the AHA dataset for writing my code because it is the simplest dataset I can work with. Simplest = significant number of baseline samples without issue of repeated measures, balanced distribution of classes, likely difference to be detected present.

p_ps = readRDS("/Users/aa370/Library/CloudStorage/Box-Box/project_davidlab/LAD_LAB_Personnel/Ammara_A/Projects/POMMS/R24_POMMS/data/processed/20221213_pomms_cleanmetadata_plant.rds")

p_ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 252 taxa and 384 samples ]
sample_data() Sample Data:       [ 384 samples by 28 sample variables ]
tax_table()   Taxonomy Table:    [ 252 taxa by 8 taxonomic ranks ]

sample_data(p_ps)

Sample Data:        [384 samples by 28 sample variables]:

c_ps = readRDS("/Users/aa370/Library/CloudStorage/Box-Box/project_davidlab/LAD_LAB_Personnel/Ammara_A/Projects/CHOICE/CHOICE_Trnl/CHOICE_20220912/ps_objects/choice_complete_20221205_clean.rds")

c_ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 571 taxa and 171 samples ]
sample_data() Sample Data:       [ 171 samples by 30 sample variables ]
tax_table()   Taxonomy Table:    [ 571 taxa by 10 taxonomic ranks ]

Preprocess

Remove problematic samples

Taking out samples with less than or equal to 0 reads and keeping only entry timepoint reduces samples from 384 down to 240.

p_ps = p_ps %>%
  subset_samples( reads > 0) %>%
  subset_samples(timepoint == "Entry")

p_ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 252 taxa and 240 samples ]
sample_data() Sample Data:       [ 240 samples by 28 sample variables ]
tax_table()   Taxonomy Table:    [ 252 taxa by 8 taxonomic ranks ]

Class distribution in this dataset: Case = 193 Control = 47

sample_data(p_ps) %>%
  as.data.frame() %>%
  as_tibble() %>%
  group_by(group) %>%
  summarise(count = n())

Merge phyloseq

combined = merge_phyloseq(c_ps, p_ps)

combined

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 722 taxa and 411 samples ]
sample_data() Sample Data:       [ 411 samples by 38 sample variables ]
tax_table()   Taxonomy Table:    [ 722 taxa by 10 taxonomic ranks ]

Aligning seq table

seqtab = combined@otu_table

seqtab_merged = dada2::collapseNoMismatch(seqtab)

combined = phyloseq(otu_table(seqtab_merged, taxa_are_rows=FALSE),
               sample_data(combined@sam_data), 
               tax_table(combined@tax_table))

combined

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 684 taxa and 411 samples ]
sample_data() Sample Data:       [ 411 samples by 38 sample variables ]
tax_table()   Taxonomy Table:    [ 684 taxa by 10 taxonomic ranks ]

#sample_data(combined)

Data transform

otu_clr = abundances(combined, "clr") %>% 
  t()

combined = phyloseq(otu_table(otu_clr, taxa_are_rows = FALSE),
                   sample_data(combined@sam_data),
                   tax_table(combined@tax_table))

rm(otu_clr)

Feature table

Extract data

Order of asv colnames in OTU table and tax table is the same.

features = combined@otu_table %>%
  as.data.frame() %>%
  as_tibble(rownames = NA) %>%
  rownames_to_column("sample")

feature_labels = c("sample")

features

Zero variance features

No features have zero variance

nzv = features %>%
  select(!sample) %>%
  nearZeroVar(saveMetrics = TRUE)

nzv %>%
  filter(nzv == TRUE)

rm(nzv)

Correlated predictors

summary(feature_cor[upper.tri(feature_cor)])

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-0.388970 -0.007175  0.037148  0.124461  0.164990  1.000000 

high_cor = findCorrelation(feature_cor, cutoff = .50) #using a more permissive cutoff because a lot of columns still left in the input table

filtered_features = features[,-high_cor]

rm(high_cor)
rm(feature_cor)

Linear dependencies

Input table

Modeling

Baseline performance

label_list = c("group", "group_shuffle")
results = loop_label(iterations =20, input, label_list, percent = 0.8, upsampling = TRUE)

[1] "will upsample using smote on training data"
[1] "group"
[[1]]
[1] "group"

[[2]]
          Reference
Prediction Case Control
   Case     626     120
   Control  134      60

[1] "will upsample using smote on training data"
[1] "group_shuffle"
[[1]]
[1] "group_shuffle"

[[2]]
          Reference
Prediction Case Control
   Case     644     159
   Control  116      21

auc_df = results[[2]]

importance_df = results[[1]]

auc_df %>%
  pivot_longer(cols = c("group", "group_shuffle"), names_to = "feature", values_to = "AUC") %>%
  ggplot()+
  geom_boxplot(aes(x=feature, y = AUC))+
  theme_classic()

t.test(auc_df$group, auc_df$group_shuffle)

rm(auc_df)

rm(importance_df)

Loop

label = "group"

importance_df = data.frame()

cummulative_external = data.frame()

auc_list = c()

current_data = input %>%
      select(-label_list[label_list != label]) %>%
      select(-sample)

for(i in 1:20) {
  # Data splitting
    train_ind = caret::createDataPartition((current_data %>%
                                              pull(label)), #split based on label column
                                     p = 0.8,
                                     list = FALSE,
                                     times = 1)
    train = current_data[train_ind,]
    test = current_data[-train_ind,]
    
    # Upsampling
    train[[label]] = as.factor(train[[label]])
    train = train %>%
        as.data.frame()
      
    formula_smote = label %>%
      paste0("~.") %>%
      as.formula()
      
    train = DMwR::SMOTE(formula_smote, data  = train)
    
    # Training model
    
    formula = label %>%
      paste0("~.") %>%
      as.formula()
    
    fit_control = caret::trainControl(method = "CV", # 5-fold cross validation
                                 number = 5,
                                 summaryFunction = prSummary,
                                 classProbs = T)
    
    model = caret::train(formula, #using formula specified 
              data = train,
                  method = "rf",
                  trControl = fit_control,
                  verbose = FALSE, 
                  metric = "AUC")
    
    # Model predictions
    
    predictions = predict(model, newdata = test, type= "prob")
    
    predictions$obs = test %>%
      pull(label) %>%
      as.factor() 
    
    # Evaluate model
    
    eval = MLeval::evalm(predictions, showplots = FALSE, silent = TRUE)
    
    auc = eval$stdres$`Group1`[13,1]
    
    auc_list = append(auc_list, auc)
    
    #Extract feature importance

    feature_importance = caret::varImp(model$finalModel) %>% # Ben calls this importances instead of final model but the rankings are exactly the same, only number value scale is different.
      arrange(desc(Overall)) %>%
      as.matrix() %>%
      t() %>%
      as.data.frame()
    
    feature_importance$auc = auc
    
    feature_importance$iteration = i
    
    feature_importance$label = label
    
    importance_df = rbind(importance_df, feature_importance)
    
    #Predicting on external dataset
    
    external_predictions =  predict(model, newdata = external_test, type= "prob")
    
    external_predictions$sample = external_test$sample
    
    external_predictions$auc = auc
    
    external_predictions$pred = predict(model, newdata = external_test)%>%
      as.factor ()
    
    external_predictions$iteration = i
    
    cummulative_external = rbind(cummulative_external, external_predictions)
    
}

auc_list

 [1] 0.61 0.61 0.48 0.76 0.74 0.32 0.82 0.46 0.68 0.70 0.78 0.68 0.70
[14] 0.70 0.56 0.57 0.42 0.68 0.65 0.57

importance_df

cummulative_external  

NA

Plotting

AUC plot and t-test

data.frame(auc = auc_list)%>%
  ggplot()+
  geom_boxplot(aes(y = auc))+
  theme_classic()

Save external predictions

predictions

write_csv(save, "/Users/aa370/Library/CloudStorage/Box-Box/project_davidlab/LAD_LAB_Personnel/Ammara_A/Projects/Machine_learning/LAD_lab_ML/code/pomms_to_choice_predictions.csv")

variable importances

write_csv(importance_df, "/Users/aa370/Library/CloudStorage/Box-Box/project_davidlab/LAD_LAB_Personnel/Ammara_A/Projects/Machine_learning/LAD_lab_ML/code/pomms_to_choice_feature_importance.csv")